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An extended Hamiltonian approach to conduct isothermal-isobaric molecular dynamics simula- 
tions with full cell flexibility is presented. The components of the metric tensor are used as the 
fictitious degrees of freedom for the cell, thus avoiding the problem of spurious cell rotations and 
artificial symmetry breaking effects present in the original Parrinello-Rahman scheme. This is com- 
plemented by the Nose-Poincare approach for isothermal sampling. The combination of these two 
approaches leads to equations of motion that are Hamiltonian in structure, and which can therefore 
be solved numerically using recently developed powerful symplectic integrators. One such inte- 
grator, the generalised leap-frog, is employed to provide a numerical algorithm for integrating the 
isothermal-isobaric equations of motion obtained. 
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I. INTRODUCTION 



Over the last few decades advances in techniques, mod- 
els rtg|d computer power have made simulation meth- 
odstrHu an indispensable aid to research in condensed 
matter, molecular physics and chemistry, and in mate- 
rials science. By emulating the conditions under which 
experiments are carried out and the interactions between 
the components of the system as closely as possible, sim- 
ulations can often provide information that is not directly 
attainable from the experiments themselves, and can help 
to interpret the empirical observations. 

There are two large groups of simulation methods 
that are capable of taking thermal effects into account, 
namely the Monte Carlo (MC) and the Molecular Dy- 
namics (MD) methods. In this work it is the latter that 
will be of interest. Conventional MD consists of numeri- 
cally integrating the classical equations of motion for an 
ensemble of atoms representative of the system of inter- 
est. Assuming ergodicity, a sufficiently long trajectory 
will sample the whole of the accessible phase space un- 
der the conditions of the simulation, and a time aver- 
age over such a trajectory of any property of the sim- 
ulated system will provide an estimate of the value of 
that property in the real system under the same con- 
ditions. The classical equations of motion conserve the 
total energy, so this procedure simulates the system of 
interest in the microcanonical (constant number of par- 
ticles, N, constant volume, V, and constant energy, E, 
or NVE) ensemble. However, experiments arc most fre- 
quently conducted under conditions of constant temper- 
ature, and sometimes also constant pressure, conditions 
which correspond to the canonical (NVT) or isothermal- 
isobaric (NPT) ensembles respectively. It would there- 
fore be desirable to have MD methods that were capable 
of sampling these ensembles also. In a seminal paper, An- 
dersenQ proposed two different methods for conducting 
MD simulations in the canonical and iso-shape isobaric- 
isoenthalpic (NPH, where H is the enthalpy) ensembles, 
methods which could be combined to perform simulations 



sampling the isothermal-isobaric ensemble. Andersen's 
NVT MD method involved a series of stochastic colli- 
sions which changed the velocity of a randomly chosen 
atom in the system to one generated from a Maxwell- 
Boltzmann distribution of velocities at the desired tem- 
perature. Such a series of collisions changes the total en- 
ergy of the system so that it fluctuates around its equilib- 
rium value as prescribed by the canonical ensemble. As 
for the constant-pressure (isobaric-isoenthalpic or NPH) 
ensemble, Andersen showed that it could be sampled by 
incorporating the volume as a new degree of freedom in 
the classical equations of motion, with a fictitious mass 
and velocity associated to it. 

Andersen's idea of extending the dynamics of a physi- 
cal system by including fictitious degrees of freedom has 
proved to be extremely powerful. Andersen himsclna pos- 
tulated that using this method it might be possible to 
perform NVT dynamics in a non-stochastic fashion. In- 
deed, such a method was put forward by NoseoQ, and 
later modified by HooverQ.r,Xhe same idea was to be used 
by Parrinello and RahmanQ'El for flexible-cell (as opposed 
to iso-shape) constant-pressure MD, and by Car and Par- 
rinellal3, who combined in an ingenious way the classical 
dynamics of ions with the fictitious dynamics of the elec- 
tronic orbitals expanded in a plane- wave basis set, within 
a Density Functional Theory description of the electronic 
structure. This was done in such a way that as the ions 
moved, the electronic orbitals followed adiabatically, so 
that the system remained in the ground state, or suffi- 
ciently close to it, allowing Car and Parrinello to carry 
out ab initio MD for the first time. 

Although the extended dynamics approach has proved 
very useful for sampling ensembles beyond the micro- 
canonical one, over time a number of shortcommings have 
been detected, both with the Nose-Hoover method for 
NVT sampling and the Parrinello-Rahman method for 
NPH simulations. The method originally proposed by 
Nose involved a time transformation, as a result of which 
the phase space is not sampled at regular time inter- 
vals. While this in itself does not pose a problem for 
the calculation of thermal averages of time-independent 
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quantities, it severely complicates the calculation of time- 
correlation functions. It is possible to|-.transform the 
equations—of motion back to real timeoQ using a non- 
canonicalEj transformation, but then the Hamiltonian 
structure is lost, i.e. the equations that result cannot be 
derived from a Hamiltonian. Although the lack of Hamil- 
tonian structure in the Nose-Hoover equations noses no 
practical difficulties (both iterativeo and explicitO inte- 
grators can be derived for these equations), it is unde- 
sirable from a formal point of view, not least because 
the usual machinery of statistical mechanics is not di- 
rectly applicable to non-Hamiltonian systems (see the 
recent work of Tuckerman and coworkersEiai-3 concern-, 
ing this point). Recently, Bond, Leimkuhler and LairdE-3 
have shown that it is possible to obtain NVT sampling 
(subject to the usual assumption of ergodicity) starting 
from Nose's original Hamiltonian, but acted upon by a 
Poincare transformation. The transformation is canoni- 
cal, and therefore preserves the form of the equations of 
motion; thus, all the recent developments in integration 
schemes for Hamiltonian systemstj (see below) can be 
used to solve numerically the resulting equations of mo- 
tion. Furthermore, the Poincare transformation is cho- 
sen in such a way that the sampling of the phase space 
takes place M regular time intervals. Thus, the scheme of 
Bond et al£3 overcomes the negative aspects associated 
with the Nose-Hoover method. 

Concerning the constant-pressure schemes, Parrinello 
and Rahman formulated their method taking as extended 
variables the Cartesian components of the simulation cell 
vectors, and constructing a fictitious kinetic energy term 
from their squared velocities . As noted by ClevelandEZl 
and by Wentzcovitchta, this choice leads to equations of 
motion that are not invariant under transformations be- 
tween equivalent cells (modular transformations) , leading 
to unphysical symmetry breaking effects. An added diffi- 
culty is that the fictitious dynamics based on this choice 
of variables results in spurious cell rotations, which com- 
plicate the analysis of the results. Numerous proposals 
to overcome, these-difficulties have appeared in the liter- 
atureOli3lijc30'E3, but here I am goins-to focus on the 
method proposed by Souza and MartinsEa, which uses as 
fictitious dynamical variables the components of the cell 
metric tensor. The metric tensor is independent of the 
cell orientation, and by phrasing the dynamics in terms 
of it, spurious cell rotations simply do not appear. Fur- 
thermore, the fictitious kinetic energy associated with the 
metric tensor can easily be constructed such that it is in- 
variant with respect to modular transformations. 

Important advances have also been achieved in the un- 
derstanding and design of integrators for dynamical sys- 
temsEJ'Ell Classical Hamiltonian dynamics is time re- 
versible (if the Hamiltonian is an even function of the 
momenta) and symplectic, i.e. it preserves the sum of 
areas spanned by the vector products dpi x dqi, the area 
element around the point (p,, q{). It is desirable that nu- 
merical integrators for Hamiltonian equations of motion 
respect the symmetries underlying Hamiltonian dynam- 



ics, such as time reversibility and symplecticness, as then 
one can be more confident that the discrete-time solution 
will resemble more closely the exact solution. Thus con- 
siderable efforts have been devoted to the development 
of integrators which comply with these requirements, and 
modern simulation techniques should take advantage of 
the progress achieved in this field. 

Constant-pressure algorithms such as the method of 
Parrinello-Rahman and its variants sample the isobaric- 
isoenthalpic ensemble. However, this ensemble is not 
very common, nor indeed is it very convenient. In many 
circumstances it is desirable to have control over the av- 
erage temperature, and it is therefore preferable to per- 
form a simulation which samples the isothermal-isobaric 
ensemble. This can be achieved by combining one of 
the constant-pressure algorithms with a Nose-type ther- 
mostat. The purpose of this paper is to describe how. 
the metric-tensor based scheme of Souza and MartinsEj 
for constant-pressure simulations can be_combincd with 
the Nose-Poincare scheme of Bond et aZilj to provide an 
extended Hamiltonian for isothermal-isobaric MD sim- 
ulations with full-cell flexibility, incorporating the ad- 
vantages of both methods over the original schemes of 
Nose-Hoover and Parrinello-Rahman. Similar, schemes 
have been recently reresented by Sergi et alts and by 
Sturgeon and LairdcJ, but only in the case of isotropic 
cell fluctuations. It is shown how the equations of mo- 
tion that result in the present scheme can be integr. 
numerically using the generalised leap-frog methodOl 
leading to a symplectic, highly stable algorithm. This 
NPT MD algorithm is then illustrated with a series of 
realistic test cases, namely NPT simulations of diamond 
and crystalline silicon. 

The structure of the paper is as follows: in section II 
the necessary background on the metric-tensor constant- 
pressure and Nose-Poincare isothermal formalisms is re- 
viewed, both schemes are then combined, and a recipe for 
integrating the resulting equations of motion numerically 
is provided. The methodology is then applied to the test 
cases in section |IH|, and the conclusions are summarized 
in section 
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II. METHODOLOGY 

Classical Lagrangian and Hamiltonian mechanicsEll 
provide elegant theoretical frameworks in which to set 
up equations of motion for the system under study, using 
the variables or coordinates most convenient in each situ- 
ation. The generalized leap-frog integration scheme that 
will be used below is most readily applied to Hamiltonian 
equations, and therefore the Hamiltonian formulation of 
classical mechanics will be used in what follows, although 
the Lagrangian treatment is of course totally equivalent. 

When considering a system in which both the volume 
and shape of the simulation cell are allowed to evolve 
in time, it is helpful to take into account the covari- 
ant/contravariant character of the dynamical variables. 
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Therefore, in the following, the covariant or contravari- 
ant character of the different variables that will be used 
shall be indicated by means of a subindex or superindex, 
respectively. 

Let us now review the set of dynamical variables which 
will be used to specify the state of the system. Firstly, 
it will be most convenient to use lattice coordinates, q l , 
which give the position of atom i relative to the simula- 
tion cell. The dynamics of atoms is then fully accounted 
for by considering the momenta p i; conjugate to q*. Lat- 
tice coordinates are related to the usual Cartesian ones, 
fj, through 



Hq', 



(1) 



where H is a 3 x 3 matrix formed by the simulation cell 
vectors a Q , (a — 1, 2, 3) in columns. Since the cell vectors 
a a are linearly independent, matrix H can be inverted, 
and its inverse, H _1 , has row a equal to b Q , the recipro- 
cal vector of & a in the following sense: 



(2) 



The calculation of interatomic distances when atom posi- 
tions are specified in terms of lattice coordinates is given 
by: 



rij = V(q l -q J )G(q l -q J ): (3) 

where G is the metric tensor, with elements 



G, 



aB 



a,, 



a 



(4) 



The volume of the simulation cell is also given by the 
metric tensor as V ce u = %/det G. The set of vectors b" 
also define a metric tensor, G™' 3 , which is reciprocal to 
the previous one. 

Souza and MartinsE3 have shown that the metric ten- 
sor constitutes a very convenient dynamical variable for 
constant-pressure MD simulations. Firstly, G is invari- 
ant under cell rotations; the orientation of the cell is ir- 
relevant, and thus spurious cell rotations do not appear 
during the dynamics. Secondly, it is easy to set up a 
fictitious kinetic energy term associated with the met- 
ric tensor (see below) which is invariant under modular 
transformationsEa (i.e. transformations between the dif- 
ferent possible cells compatible with the periodicity of 
the system) . This avoids artificial symmetry breaking ef- 
fects. Thus, spurious cell rotations and symmetry break- 
ing effects, which appeared in the original Parrinello- 
RahmantM constant-pressure algorithm, are naturally 
avoided in this formalism. Following Souza and Mar- 
tins, each metric tensor component G a p has a conjugate 
momentum P™' 3 , and the fictitious kinetic energy term 
associated to the dynamics of the metric tensor is 



K a = 



pa pB 

B a 

2M G det G ' 



(5) 



where the sum over repeated indices is implied. Here 
Mo is a fictitious mass, but the total effective mass is 



Mq det G, which varies with the cell volume. While it 
would be possible to use a constant fictitious mass, this 
form has the particularity of reducing the kinetic energy 
expression in Eq. (|^) to the same form as in Ander- 
sen's constant-pressure scheme in the case of iso-shape 
cell fluctuations. 

In the case of hydrostatic external pressure, the po- 
tential energy term associated to the metric dynamics is 
simply 



U G = Text Vcell = V ex tVdetG. 



(6) 



Souza and Martins went on to show that the case of an 
anisotropic external stress can also be contemplated, if a 
potential energy term of the form 



Uc 



is included, where <rf " t are the components of the external 
stress in contravariant lattice coordinates. 

The combined dynamics of atoms and metric tensor de- 
scribed so far conserves the enthalpy (the generalized en- 
thalpy of Thurstoncj, in the case of constant anisotropic 
external stress), and samples the isobaric-isoenthalpic 
(NPH) ensemble. It is therefore desirable to combine 
the dynamics of the extended system of atoms and met- 
ric tensor with a device that enables sampling of the 
isobaric-isothermal (NPT) ensemble. Souza and Martins 
used stochastic Langevin dynamics in the examples re- 
ported in their work. A different strategy will be pursued 
here, which consists of coupling the dynamics of atoms 
and metric tensor with a Nose thermostat, as described 
in what follows. 

Canonical (NVT) MD simulations have been usu- 
ally undertaken by means of the so called Nose-Hoover 
method. Recently, however, Bond et al.lB have provided 
an alternative scheme, which also samples the NVT en- 
semble, but has the additional advantage of being Hamil- 
tonian in structure, thereby permitting the use of sym- 
plectic numerical integrators. This is achieved by per- 
forming a Poincare transformation on the original Nose 
Hamiltonian, i^Nosc, which results in 



1 



./3a 



G a 







(7) 



osc-Foincarc 



H ), 



(8) 



where Hq is a suitably chosen constant, and i?Nosc is 
given by (in Cartesian coordinates) 



H^osc = 2_j 



2rrn S 2 



U{v) 



p2 

2M S 



9 k B T ext In S. (9) 



Here S is the position variable of the thermostat, a 
strictly positive quantity, Ps its conjugate momentum, 
g is the number of degrees of freedom of the physical 
system (i.e. excluding extended or fictitious dynami- 
cal variables), k B is Boltzmann's constant, and T ex f is 
the temperature of the thermostat. Bond and coworkers 
have demonstrated that, under the assumption of ergod- 
icity, the Nose-Poincare Hamiltonian generates dynamics 
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that sample the canonical ensemble, as desired. Recently, 
Sturgeon and LairoO have extended the Nose-Poincare 
Hamiltonian with an Andersen barostat, which imple- 
ments iso-shape cell fluctuations, thereby sampling the 
isothermal-isobaric (NPT) ensemble appropriate for non- 
crystalline systems. For crystalline solids, however, it is 
necessary to allow fluctuations of the cell shape as well as 
of its volume, which can be done by means of the metric 
tensor dynamics described above. The scheme presented 
here can therefore be regarded as a generalization of the 
method reported by Sturgeon and Laird to full-cell dy- 
namics NPT MD simulations. 

By combining the Np&e-Poincare Hamiltonian [Eq. (^)] 
of Bond and coworkerst-3 with the metric-tensor constant- 
pressure scheme of Souza and Martina^, one arrives at 
the following Hamiltonian: 



In this equation, lower case labels p and q refer to atomic 
momenta and position variables of the physical system 
of interest, while labels in upper case refer to fictitious 
degrees of freedom (thermostat or barostat). Latin in- 
dices label atoms, always appear first and are always 
used as subindices, while Greek indices label compo- 
nents of tensors, and if used as subindices appear after 
the atom label. Sums over atoms, as in the atomic ki- 
netic energy term, are written out explicitly; otherwise 
the summation over repeated indices of tensorial quan- 
tities is implied. Thus, Pi a is the covariant a compo- 
nent of the momentum of atom i, while P a o is the mixed 
(contravariant-covariant) component of the second-rank 
tensor P a p formed by the momenta associated to the 
components of the metric tensor G. The constant Hq 
is to be chosen so that £/npt has a value of zero. 



H 



NPT 



5' 



E 



PtaPT 

2 m, S 2 



"W(q, G) 



pa p/3 
r g r a 

2M G det G 
1 



V ext VdetG + tofc G aP +(10) 



P 2 
2Mc 



gk B T ext InS - H 



Using the standard rules of Classical MechanicsEl, and 
with the help of the relations d det G /dG a p — G^ a det G 
and dG Xf */dG af 3 = -G Aq G^, it is straight forward to 
obtain the following equations of motion: 



Qi 

Pia 
G a p 
pa/3 



s 

Ps 



S 



-s 



= s- 



dU 



Pap 



M G det G 



= -S 



dU 



Pi Pi 



(11a) 

(lib) 

(11c) 



P px G Am P^ a 



dG aP ^2m,5 2 Mg det G 



= S 



Ps 



M s 

\ ^ Pia Pf 

/ — J rrii S 2 



g ksTe.. 



-V ext VdetG 



AH 



2 M G det G 



Q0a 



1 0a 

2 ^ext 



(lid) 

(lie) 
(llf) 



where the dot indicates a time derivative, and 



AH = 



sr^ p^ Pi 



^ 2 rm S 2 



■W(q,G) 



2 M G det G 



+V ext VfetG+-<T^ t G af3 



gk B T ext In 5 - H . 



(12) 



Let us now turn to the question of how to obtain a 
numerical scheme to integrate these equations of motion. 
Of the different schemes that are possible, I have cho- 
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sen to use the generalized leap-frog scheme (GLF)@ 
which is simple, and since the system is Hamiltonianlij, 
the numerical procedure that results from applying-tb 



GLF to Eqs. (Ill) is symplectic and time-reversibL 



Let us briefly recall how the GLF works: let Q, P repre- 
sent posit ion and momentum variables respectively. The 
Eqs. (Hi) can be generally put in the form 



P 



G(P,Q), 
--F(P,Q). 



(13) 
(14) 
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The GLF scheme consists of first propagating the mo- 
mentum variables half a time-step forward from the ini- 
tial point in phase space: 

P(t + \/2 5t) = P{t) + 

St F[P(t+ 1/2 St), Q(t)]/2, (15) 

where St is the time step. Note that the new momenta ap- 
pear on both sides of the equation, which can lead (and 
in general, does lead) to implicit equations for the up- 
dated momentum variables. Next, the position variables 
are updated a full time-step: 

Q(t + 5t) = Q(t) + 

St {G[P(t + 1/2 St), Q{t)}+ (16) 
G[P(t+ 1/2 6t),Q(t + 6t)]}/2. 



Again, due to the fact that Q(t + St) appears on both 
sides, in the general case one obtains an implicit equation, 
like for the momentum variables at half step. Once Q(t + 
St) has been obtained, new forces can be calculated, and 
thus it is possible to bring up the momenta to full step: 

P(t + St) = P(t+ 1/2 Si) + 

St F[P(t + 1/2 St), Q{t + St)]/2. (17) 

The scheme can then be iterated by setting t = t + St 
and returning to Eq. (|l5|). 

Applying the GLF to the NPT equations of mo- 
tion (|llj ) leads to the following numerical scheme: 



Pia,l/2 



1/2 



P.S,l/2 

Si 

G a B,l 



P 



a/3 



PietA 



P 



a0 



\stS Q 



dU 



8G, 



a/3,0 



E 



Pi, 1/2 PiS/2 

2m t Sl 



p/3A ^ pAia 
^1/2 ^A^.O 



J72 

Mq det G 



/2 



Ba 



Ps,o + 



E 



Pux,l/2 Pi,l/2 



i k B T, 



BJ-ext 



A^(qo,G ,5o, Pl /2i Pg,1/2i Ps,1/2 



Si 



S,l/2 



1 / Gp\.Q P]_?2 Gpafi 
G a Bfi + ^St S — 

2 V Mq det G 



Si 



Mq detGi 



Co 



i,l/2 



2 \m, b rriibi 



Ps,l = 



1 



P S ,1/2 + ^St 



E 



Pia.l/2 P iA/2 



gkBT e .. 



Aff (qi, Gi, Si, Pi/2, P G ,l/2, Ps,x/i)] 



P i/2 - 2 St Sl 



1 



OU ^ Pi,l/2Pi,l/2 r i/2 ^1/2 



E : 



0<?a/9,i ^ 2mi5? ' M G detGi 

'A pP \ 



pA p/J 

M/2 J A, 1/2 ) £,,3 a 



1 /3a 



2 «9g" 



(18a) 



(18b) 

(18c) 
(18d) 

(18e) 

(18f) 

(18g) 

(18h) 
(18i) 



In these equations, subindices of 0, 1/2 and 1 indi- be taken at zero, half or full tim e step from the initial 
cate that the corresponding dynamical variable should point in phase space. Eqs. (|l8i|a-c) propagate the mo- 
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menta forward in time by half a time step; Eqs. ( |l8ip -f) 
advance the position variab les the whole length of the 
time step, and finally Eqs. ( 18i g-i) complete the updat- 
ing of the momenta to full time-step. The order in which 
Eqs. (18i) arc written corresponds to the order in which 
they must be implemented in a computer code. Note 
that, as indicated above, som e of the equations are im- 
plicit. In particular, Eq. ( |l8i| b) is implicit in the metric 
tensor momenta at half time step, and must be solved 
iteratively. This is also the case for the metric tensor 
components themselves, in Eq. ( |18i| e). Two strategies 
can be adopted for solving these equations: a Newton- 
Raphsor£3 procedure can be easily applied, but a simple 
iterative scheme seems to work equally well and is even 
simpler to program. In the example applications reported 
below, the iterative scheme was used. This consisted of 
using as initial guess for the half-step metric-tensor mo- 

1/2J 



menta in Eq. (18ib) (Pf%) their values at the start of the 



molecular dynamics time step (i.e. Pq )■ The resulting 



values for Py 2 are then fed again into Eq. (18ib), thus 
obtaining a new estimate for the correct values at half- 
step. This procedure is iterated until the absolute values 
of the differences of Py 2 found in two successive itera- 
tions of the procedure differ by less than 10~ 7 . Since, for 
a sufficiently small time step, the values of Pq 13 are close 
to those of P"jl, this procedure converges very rapidly, 
requiring only a few iterations. Exactly the same pro- 
cedure was followed for solving Eq. (|l8je). It s hou ld 
be emphasized that the implicit nature of Eqs. (18ib) 



and (|18je) does not have any significant impact on per- 
formance; the iterative procedure is very fast, and does 
not require re-evaluation of the energy, force or stress; 
however, because of the fact that a convergence criterion 
must be used, the exact time-reversibility of the equa- 
tions of mot ion [Eqs. (Hi)] is lost in the numerical scheme 
of Eqs. (18i). Nevertheless, the use of a strict convergence 
criterion makes the sche me ti me-reversible within numer- 
ical accuracy. Equation (18ic) also deserves attention; it 



is quadratic in Pg.i/2i an d care must be taken to choc. 



the right root and avoid numerical cancellation errors^ 



Equations (18i) are quite simple to incorporate into 
existing MD codes. The only additional information re- 
quired from what is conventionally computed in an MD 
program are the derivatives of the potential energy with 
respect to the metric tensor components, ■ These 
can be obtained taking into account that inter atomic 
distances depend on the components of the metric ten- 
sor, as indicated in Eq. (||). However, if the program al- 



ready computes the derivatives 



au 



where are the 



components of the Cartesian strain tensor, this informa- 
tion can be converted into the required derivatives in a 
straight forward fashion, as the strain and metric tensors 
are related througrU 



where Hq is the inverse of the undistorted cell matrix, 
and 1 is the unit matrix. Then, by simple application of 
the chain rule, one obtains the following relation: 



dU 



[(Ho'fGH^-l], 



(19) 



dG afi -V H » 1)aX ^ H » 1)T ^ (20) 
where again summation over repeated indices is implied. 

III. APPLICATIONS 

The metric tensor Nose-Poincare NPT method de- 
scribed in section II has been implemented in a computer 
program which performs MD-simulations of systems de- 
scribed with a Tight Binding! 3 ^ (TB) total energy model. 
In this section, after a brief description of the model used, 
I will report the results of tests of the methodology and 
some examples of applications directed at demonstrating 
its usefulness. 

A. Model 

The NPT algorithm described in this paper applies 
equally well to any atomistic model from which forces 
and stresses, as well as the total energy, can be obtained, 
and therefore any such models can be used in conjunction 
with it. Atomistic models can be broadly classified into 
three groups, namely empirical ■potential methods, semi- 
empirical total energy methods and first principles or ab 
initio total energy methods. Each group of methods has 
its strengths and weaknesses; for example, empirical po- 
tential methods are computationally cheap, but ad hoc 
in nature, and therefore of limited reliability and trans- 
ferability. First principles methods, on the other hand, 
rest on firm theoretical grounds, but are orders of mag- 
nitude more demanding computationally than empirical 
potentials. Added difficulties here are the bad scaling 
of the computational cost with the size of the system 
under study (usually, for electronic structure methods 
the cost scales as iV 3 or worse, where N is the number 
of atoms, altho ugh c onsiderable progress towards linear- 
scaling methodsEilH has been achieved in recent years), 
and, of specific relevance to constant-pressure simula- 
tions, the accuracy of the total energy and its deriva- 
tives depends on the volume of the simulation cell, unless 
highly converged integrationS|-wer the irreducible cell of 
the reciprocal lattice are usedl23. 

In the applications that follow, I have chosen to use an 
approximate Tight Binding (TB) modelEJ, which corre- 
sponds to the group of semi-empirical methods discussed 
above. This choice is motivated by several considera- 
tions. Being in between the extremes of empirical poten- 
tials and first principles methods, it shares some of the 
advantages (but also some of the disadvantages) of both. 
Its computational demands are very modest; it is based 
on a quantum treatment of the valence electrons, even 
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FIG. 1: Contributions to the conserved quantity Hnpt 
[Eq. ([n])] for a carbon system in the diamond structure (54- 
atom cell) as a function of time. The curve labeled physical 
energy is the energy of the atomic sub-system (kinetic plus 
potential energy). 



if at a rather simplified level, and therefore it incorpo- 
rates the essential features of the quantum nature of the 
chemical bond. In spite of their simplicity, there exist 
TB models which are capable of surprising accuracy in 
their predictions. In particular, the applications reported 
below have been carried., out using a TB model due to 
Porezag and coworkerso It goes beyond conventional 
TB models in that it incorporates the non-orthogonality 
of the basis set, which is usually assumed to be orthogo- 
nal, is constructed on the basis of Density Functional cal- 
culations employing the same basis set, and the range of 
the hopping integrals extends beyond the nearest neigh- 
bor distance, which is the range used in most TB models. 
Additional details on-this model can be found in the pa- 
per by Porezag et alS3. 



B. Diamond 

To illustrate the stability an d ac curacy of the integra- 
tion scheme embodied in Eqs. ( |18i| a-i), a simulation of di- 
amond at GPa external pressure and a temperature of 
1000 K was performed. This is well below the Debye tem- 
perature of diamond (2340 K), and therefore it is strictly 
speaking not justified to perform a classical MD simula- 
tion at this temperature. Nevertheless the aim here is to 
test the methodology, and not to extract any conclusions 
on the physics of diamond at 1000 K. The starting con- 
figuration of the system consisted of 54 atoms at their 
equilibrium positions in a cell of edge length 14.33 bohr 
(7.58 A), with the edges forming 60° angles. Initial veloc- 
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ities were chosen randomly from the Maxwell-Boltzmann 
distribution at the desired temperature, and modified to 
eliminate any translation of the center of mass. The ficti- 
tious mass of the Nose thermostat was chosen to be equal 
to the mass of a carbon atom, while for the barostat 
a mass of 10 au was used. Both thermostat and baro- 
stat were assigned zero initial momenta, and the value of 
Hnpt at zero time fixes the value of Hq. 

The total energy of the system [W(q, G) in Eq. jll])] 
was calculated using the TB model of Porezag and 
coworkersEil described above. The generalized eigenvalue 
equation that results for each system configuration ac- 
cording to this model was solved using direct diagonal- 
isation, and a set of 4 reciprocal lattice vectors chosen 
according to the Monkhorst-Pack schema^ were used, 
which were sufficient to converge the total energy better 
than 0.04 meV/atom. The length of the time step was 
set to 1 fs, and the simulation was run for a total of 10 ps. 

In Fig . (0) the different terms contributing to Hnpt 
[Eq. ([n])] have been plotted as a function of time for 
the first 0.5 ps of the simulation. The physical energy is 
the sum of the kinetic and potential energy terms, cal- 
culated from the atomic momenta and from U{(\ 1 G), re- 
spectively. The thermostat and barostat contributions 
include both the kinetic and potential energy terms as- 
sociated to each of these fictitious degrees of freedom 
(the barostat energy has been scaled by a factor of 100 
so that it can be appreciated on the graph). Note that, 
as expected, the total energy of the atoms (their kinetic 
plus potential energy, the physical energy), is not con- 
served, contrary to what would happen in a conventional 
microcanonical (NVE) MD simulation. In the present 
case, it is -/Znpt that is the conserved quantity, and as 
can be seen from Fig. (jlj), this is indeed ap prox imately 
conserved by the numerical scheme of Eqs. (18ia-i). To 
judge how well i?NPT is conserved during the dynamics, 
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FIG. 3: Cell vector lengths as a function of time for the 54- 
atom diamond cell. 



FIG. 4: Cell angles as a function of time for the 54-atom 
diamond cell. 



it has been plotted in Fig. (g) for the whole length of the 
simulation. Except at the very first stages of the calcu- 
lation, the deviations of Hnpt from its mean value are 
smaller than 0.0005 hartree, and its standard deviation 
has been computed to be 0.00015 hartree. No drift is ob- 
served during the trajectory; simulations of up to 50 ps 
(see below) were also carried out, and again, no appre- 
ciable drift was observed, testifying to the stability of the 
method. 

Figs. (||) and (f|) illustrate the time evolution of the 
super-cell edge lengths and angles during the first 5 ps 
of the simulation, respectively. The edges start at time 
zero having equal lengths, corresponding to that of a 54- 
atom cell at K, but as the simulation proceeds each 
cell parameter evolves separately. Note how the average 
values of the moduli of the cell vectors settle at a higher 
value than the K one, namely at 14.41 bohr (obtained 
by averaging over the whole length of the simulation). 
This is due to the thermal expansion of diamond from 
to 1000 K, though it should be emphasized that one 
must not expect an accurate estimation of the thermal 
expansion in this case, partly because of the complete 
neglect of quantum effects, which, as indicated above, 
are important below the Debye temperature. The cell 
angles also evolve independently shortly after the start 
of the simulation, but contrary to what happens with 
the cell vector moduli, which evolve to a different mean 
value, the cell angles oscillate around their initial value 
of 60°, with their instantaneous values remaining within 
±0.4°, indicating that, although the cell expands due to 
the thermal motion of the atoms at 1000 K, it does not 
change its shape. 

Fig. (||) illustrates the time evolution of the cell vol- 
ume and the internal pressure during the first 5 ps of 
the simulation. As expected, these two magnitudes dis- 



play opposite behavior, in the sense that when the vol- 
ume is lowest, the internal pressure is highest, and vice 
versa. Like the cell vector moduli [Fig. (H)], the volume 
expands from the K value to a slightly larger value ap- 
propriate to the temperature of the simulation, around 
which its value oscillates. The internal pressure, on the 
other hand, oscillates around GPa, the value fixed for 
the external pressure in this simulation. At the end of 
the trajectory, the thermal averages of the internal pres- 
sure and the temperature of the system were computed to 
be 0.2 GPa and 999.7 K respectively, in good agreement 
with the imposed external values. 
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FIG. 5: Cell volume (left ordinate axis) and internal pressure 
(right ordinate axis) for the 54-atom diamond cell. 
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FIG. 6: Silicon cell parameter as a function of temperature 
obtained from NPT MD simulations at GPa external pres- 
sure and a series of temperatures. 



C. Thermal expansion of silicon 

Silicon has a smaller Debye temperature than dia- 
mond, c.a. 640 K, and I will therefore use classical MD 
to study its thermal expansion above this temperature. 
A series of NPT simulations at GPa external pressure 
and different temperatures have been carried out. The 
total length of the simulation was 50 ps, using a time 
step of 1 fs. The external temperature of the simulation 
was varied in steps of 50 K between 700 and 1600 K, and 
at each such temperatures a simulation was conducted. 
In all of them the value of -Hnpt was accurately con- 
served, and no appreciable drift in its value was observed. 
The variation of the silicon cell parameter with tempera- 
ture is illustrated in Fig. (^J). Qualitatively, two approx- 
imately linear behaviors can be observed, between 700 
and 1150 K, and from 1150 to 1600 K, with the second 
range of temperatures having a slightly lower slope. As 
can be seen, the behavior is not very smooth, indicating 
that there is still a degree of statistical noise in the ther- 
mal averages of the cell parameter computed from these 
simulations. I will dwell on the possible causes of this 
below. 

In principle, from the data shown in Fig. (^|) it would 
be possible to calculate the thermal expansion coefficient, 
a, of silicon in the range of temperatures considered, but 
the statistical inaccuracies present in the data make this 
a difficult task. The experimental value at 1000 K is 
Qex-p = 4.3 x 1Q _6 K , while from the data in the figure 
one can estimate a value of a ca i c = 8.1 x 10~ 6 K _1 at 
this temperature. The thermal expansion coefficient is 
extremely sensitive, and given the statistical uncertain- 
ties present in the results, the calculated value can only 
be considered a rather crude estimate of the value pre- 



dicted by the TB model used; a more accurate estimate 
could be closer to the experimental value. 

It is instructive to consider the possible causes of the 
poor statistics observed in Fig. (^). It is certainly not 
due to inadequate conservation of i?NPT, which is suffi- 
ciently well conserved throughout all the simulations. A 
more likely explanation is that the dynamics generated 
is not sampling the NPT ensemble with sufficient effi- 
ciency. Inefficient sampling can occur because some de- 
grees of freedom do not easily exchange energy with the 
rest of the system, a situation which takes place when 
there are largely different frequencies present. Fig. (||) 
lends some weight to this consideration. There it can 
be seen that the volume is oscillating quasi-harmonically 
with a single dominant frequency. The internal pressure, 
however, has two dominant frequencies: a high one, re- 
flecting the thermal vibrations of the atoms, and a lower 
one, with the same frequency as the volume. Sampling 
efficiency could be increased by reducing the difference 
between the high frequency oscillations in the internal 
pressure and the frequency of the volume motion, which 
can be achieved using a lower barostat fictitious mass. It 
should be pointed out that, formally at least, the NPT 
ensemble is sampled independently of the values used for 
the masses of the fictitious degrees of freedom, provided 
the dynamics is ergodic. The sampling efficiency, how- 
ever, does depend on the values chosen, and therefore 
this choice must be made with care. 



D. Silicon under uniaxial external stress 



In this final example, the capability of the method to 
cope with the simulation of systems subjected to non- 
hydrostatic external pressures will be illustrated. Again, 
a 64 atom silicon supercell in the diamond structure was 
used as test case. The external temperature was fixed 
at 1000 K, and a series of simulations were performed 
applying an external stress in the [100] direction, varying 
from -5 to 5 GPa in steps of 1 GPa. In the sign convention 
used here a positive sign indicates a compressive pressure. 

In Fig. (0) the variation of the cell axis as a function of 
the applied stress is shown. The supercell length in the 
direction of the applied stress has a negative slope, and it 
starts at an expanded value, larger than the equilibrium 
value at the same temperature and zero stress, reaching 
a value below this when the applied stress is compres- 
sive. The other two cell lengths have positive, though 
smaller, slopes, expanding as the first dimension is com- 
pressed. Both |b| and |c| show nearly identical behavior, 
as expected, given the symmetry of the system. This fig- 
ure illustrates how the elastic constants of materials (in 
this case the Poisson ratio) could be evaluated at finite 
temperatures using this methodology. 
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IV. CONCLUSIONS 

In this paper I have presented a stable, symplcctic 
algorithm for integrating the Hamiltonian equations of 
motion pysulting from the combination of the Souza and 
MartinsEa metric tensor-based constant-pressure scheme 
and the Nase-Poincare thermostat scheme of Bond and 
coworkersEj. The dynamics generated by these equa- 
tions samples the isothermal-isobaric (NPT) ensemble 
with full-cell flexibility. Conditions of non-hydrostatic 
external pressure can also be simulated. The numeri- 
cal scheme advocated here is easy to implement in exist- 
ing molecular dynamics codes. The capabilities of this 
methodology have been illustrated with a series of nu- 
merical experiments in carbon and silicon in the diamond 
structure, using a tight binding model. 



FIG. 7: Silicon cell vector lengths as a function of the applied , , . , 

uniaxial external stress in the [100] direction from simulations 
at 1000 K external temperature. 
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